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ABSTRACT 



The interaction of land 5 kev copper atom (primaries) impacting the 
Oil surface of copper crystals was studied using an "n-body type" 
computer simulation program. Each primary energy was studied using 
both the Bohr and Gibson II copper on copper potentials. Primary pene- 
tration and lateral drift or spread were investigated as a function of 
impact point. 

The n-body program has greater penetration values than the binary 

13 

model used by Robinson and Oen . Drift or spread (defined as wander 
if energy loss rate is low) occurs in the impact region of transition 
from hard interaction (near a lattice atom) to soft interaction (open or 
channel regions). Penetration and spread contours are presented. 

The guidance and assistance given by Associate Professor Don E. 
Harrison,, Ir. of the U. S. Naval Postgraduate School in this study is 
gri'Vj t fully a aknowledged . 
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1. INTRODUCTION. 



The subject of crystal bombardment by energetic ions or atoms may 
be discussed from either external effects or internal interactions , 

Neither phenomenon is well understood, nor have adequate theories been 
devised to explain these interactions. One external effect is described 
by the sputtering ratio, which is defined as the ratio of ejected atoms 
from a crystal per impact atom. A more general term for crystal modi- 
fication caused by incoming energetic atoms is radiation damage. 

When incoming atoms or ions (called primaries) penetrate the crystal, 
a series of collisions may take place in which the primary energy is 
distributed to the crystal atoms. Often crystal atoms leave their sites 
and become moving atoms called secondaries. The orientation of the 
crystal, relative to the primary beam direction, is a critical parameter 
in any description of the interactions that take place. Deeper penetrations 
occur when primaries enter more open or transparent channels in crys- 
tals . 

The primary energy effectively categorizes the type of interaction 
which will occur. Heavy incoming atoms with relatively low energies 
usually interact with the atoms of the crystal through elastic collisions. 
Such atoms deviate markedly from straight line trajectories. High energy 
primaries interact by ionization and electron excitation of the crystal 
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and, as a result, lose little energy. These atoms usually follow 

1 * 

relatively straight line trajectories. 

This study is concerned with the ranges and spread of lower ener- 
gy primaries as they dissipate energy through internal elastic inter- 
actions in the crystals. Range is defined as the depth of penetra- 
tion of the primary from the surface. For purposes of this study total 
path distance and vector range from impact point to a specified point 
are not needed. Spread is defined as the radial displacement of the 
primary from the axis of its original impact line. 

These introductory sections provide a background in the recent 
experimental and machine calculation studies of low kev atoms and 
ions in crystals. The topics of sputtering and basic radiation dam- 
age are not primary objectives of this thesis. However, several 

good sources of information are recommended: a basic text for intro- 

2 

ductory study, by Dienes and Vineyard, and more advance articles 
by Wehner, et. al. Nelson, et. al. , and Harrison, et. al.^. 

Although recent experimental range studies have been performed 
with heavy ions bombarding single crystal targets; the nature of these 
experiments and the lack of knowledge concerning potential functions, 
make the comparison between laboratory and simulation data diffucult. 
Until recently, the only theoretical potential functions available were 
* All footnotes refer to the Bibliography 
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those from homonuclear atomic interaction. Thus, the computer pro- 
grams have used primaries and crystals of the same type, while in 
laboratory experiments the primaries and crystals have been differ- 
ent substances. 

Virtual leaders in the field of experimental range studies are J. A. 

Davies and his co-workers at Chalk River, Ontario. Their technique 

is both sensitive and reliable, but requires that radioactive gas ions 

be used as primaries. The approach entails: 

. . .bombarding an optically flat aluminum target with a mono- 
energetic beam of radioactive ions, and then measuring the 
depth of penetration of the incident beam by dissolving succes- 
ive uniform layers of aluminum from the target surface and measur- 
ing the amount of radioactivity in each layer. The ion bombard- 
ments were carried out using an electrostatic accelerator spec- 
ially designed for the purpose. The technique. . . (dissolving ex- 
tremely thin layers of known thickness). . .consists essentially 
of two steps :electrochemical oxidation at constant voltage in 
aqueous ammonium citrate, followed by chemical removal of the 
anodic oxide film in a phosphoric acid chromium trioxide solution. 

Due to the highly protective nature of the anodic oxide film, this 
process permits highly reproducible surface layers of metal as 
thin as to be removed. . . .The process may be repeated as of- 
ten as desired, so that it is possible to obtain the complete range- 
distribution curve from a single bombardment. ^ 

Some of the earlier polycrystaline range studies performed by Davies, 
et. al. included .7 to 60 kev Na^ ions in aluminum,^ and 2 to 600 kev 

pc 3 

Kr°° ions in aluminum and tungsten. All of the range distribution curves 
exhibited "tails" of deep penetration. This phenomenon was attributed to 
channeling of ions in preferred directions in the crystal. Davies, et. al. 
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stated, 



The magnitude. . . (of crystal lattice effects). . .particularly 
emphasizes the need for using either single crystals or 
amorphous materials in future range experiments . 8 

q 

In a more recent study the same group studied single crystals of 
tungsten using 1 to 20 kevXe 125 ions . Tungsten has a body center- 
ed cubic (b.c.c.) structure in contrast to aluminum and Cu which have 
face centered cubic (f.c.c.) crystals. General results of the single 
crystal studies were? 

. . .crystallographic effects are even larger than those found 
previously in Al; furthermore the results are in qualitative 
agreement with the theoretical prediction for b.c.c. lattice, 
i.e. the (100) and (111) are the most favoured directions for 
tunneling and the (110) and (112) are less favoured. This 
contrasts with the f.c.c. structure, where both theory and 
experiment find (110) the most favoured and (111) one of the 
less favoured. ^ 

Lutz and Sizmann^ bombarded copper (f.c.c) single crystals with 
Kr ions of 10 to 140 kev. energy. Their results agreed in general with 
Davies' single crystal data, since they showed that 110 channel shots 
had deepest penetration. Further, the range distribution curves showed 

o 

the characteristic tail of deep penetration (up to 4000A) for 110, 111, 
and 100 channels. Their 25 kev and 140 kev range distribution curves 
(as well as some of Davies' curves) are shown in figure 1. The experi- 
mental procedure used by Lutz, et. al.^ employed the removal of thin 
layers of bombarded crystal by controlled sputtering with non-radio- 
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active krypton ions. The radioactivity was then measured. 

In conjunction with the monocrystal experiments described above 
and computer simulation studies to be discussed, Lehmann and Leib- 
fried 11 have studied the behavior of primaries in channels analytical- 
ly. Their work used copper primaries in copper, as have the majori- 
ty of computer programs. Their results, oriented toward long range 
penetrations, in part were: 

Two potentials are used: an exponentially screened Coulomb 
potential after Bohr, used also in the machine calculation and 
thought to give an adequate description for relatively high en- 
ergies and small interatomic distances; and a purely exponen- 
tial potential after Born-Mayer, better suited for relatively low 
energies and large atomic distances. The maximum ranges are 
very large, for 10 kev in the order of 10^ lattice parameters for 
the Born-Mayer potential and up to 10 7 for the Bohr potential. 
Presumably, the Born-Mayer potential is a better description for 
these events. The investigation is confined to motions near 
the channel axis . 
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2. COMPUTER SIMULATION STUDIES 



This project studies the interactions of energetic atoms in crys- 
tals through the use of a digital computer simulation programoThis 
section will present a general orientation of the use of computer sim- 
ulation programs and will discuss the current work in this area. (A 
more detailed discussion of the program and model used in this study 
will be presented in section four and in Appendix II.) 

Simulation programs for any system have the following three es- 
sential ingredients. (1) The system or phenomenon has a status. 

This status is described in terms of attributes. In the case of a crys- 
tal, such attributes as position, velocity, mass and energy are essen- 
tial. (2) The status is changed or modified by an event. Events cause 
the atoms to leave their sites, gain or lose energy, etc. (3) The 
third essential element of a simulation program is time. The status and 
events are strongly intertwined in a time continuum. The control of the 
occurance of events is critical. 

Three basic approaches to simulation programs of atom-crystal in- 
teractions that have been used to date are: (l)rartiom target atoms (i.e. 
lattice effects are neglected), (2) binary collision in a lattice, and 
(3) many body collision in a lattice. The binary collision assumes in- 
teraction with one atom at a time, while the many body approach (called 

n-body) simultaneously accounts for forces from several crystal atoms. 

12 

One of the earliest computer studies was made by Gibson, et. al. 
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Their model placed the copper atoms on sites in the f.c.c. crystal 
and used an n-body approach. To simulate an infinite crystal, add- 
itional forces were added to the surface atoms. Because the poten- 
tial functions and their corresponding force functions were not 
known precisely, Gibson, et.al. used "a simple central difference 
procedure. . .which gives reasonable accuracy. Their techniques 
used f=ma in an iteration procedure. This study showed evidence 
of chains in the 100 and 110 directions. A chain is defined as the 
propagation of energy along close packed rows of atoms in a crys- 
tal without mass transport. Interatomic potential, surface effects, 
lattice defects, and other low energy (<400 ev) phenomena were 
studied. 

The random lattice mentioned above was used in an earlier study 
by Oen, et. al.* for copper primaries in solid copper. Their pro- 
gram assumed independent binary collision (accomplished by an im- 
pact parameter restriction of one half the nearest neighbors distance), 
classical scattering, and the Bohr hard sphere approximation to the 
Bohr potential. They concluded : 

It is found that neither the hard sphere approximation nor the 
inverse r squared approximation to the Bohr potential is partic- 
ularly good. To obtain correspondence with experimental re- 
sults it is found that the Bohr screening length must be in- 
creased as the atomic number of the interacting atoms increas- 
es. 1 
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A more recent program by Robinson and Oen 1 ^ (hereafter referred 
to as RO) used the binary approach again, but replaced the random 
distribution with a copper f.c.c. lattice . This program studied sev- 
eral potentials 3 the Bohr, eroded Born-Mayer, and truncated Born- 
Mayer. The Born-Mayer potential selected by RO had the same para- 
meters as the Gibson II potential, but incorporated four different 
truncation values. 

Significant among their conclusions were the following: 

(1) Ranges of primaries were strongly dependent on direction, 
and in f.c.c. crystals the order of range was (01 1)> (001)> (1 11)« 
isotropic. 

(2) Larger range values were due to many glancing collisions that 
confined the primaries to regions of low potential. This phenomenon of 
deep penetration is called channeling . 

(3) The combination of Bohr potential and 1000 ev primaries im- 
pacting in the (Oil) direction produced several primaries that moved 
from channel to channel. That is. RO found that the primaries showed 
large spread as they penetrated the crystal. Such spread was attributed 
to the nature of the Bohr potential, feee figure 2, which is a projection 
of their (Oil) channel trajectories onto the (Oil) plane. Note that three 
of the four tracks show large amounts of spread.) 
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A more recent n-body model was used by Gay and Harrison 14 
(hereafter referred to as GH) . They used an iteration procedure, 
but one that has important differences from Gibson's approach 
(see Section 4). Further, the crystal size was smaller than Gibson's, 
and the program was simplified to reduce computer time. A compari- 
son between binary collision characteristics and n-body collision 
approximations was performed. In general terms, the report showed 
that binary collision approximations approach true n-body collisions 
for small impact parameters. But at larger impact parameters the 
binary collisions have less energy transfer from primary to target, re- 
duce the primary scattering angle and slightly increase the recoil 
angle. Their work provides specific quantitative results for interac- 
tions in terms of energy exchanged, scattering angles, and recoil 
angles. 

A modified form of the GH model was used in this study. 
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3 STUDY OBJECTIVES 

This study of energetic atoms in crystals was performed using an 
n-body model approach, incorporated in a computer simulation pro- 
gram. Specifically the study objectives were? 

(1) The verification and more explicit delineation of the spread 

13 

phenomenon found by RO, (figure 2). As mentioned earlier, three 
of the four copper primaries shown in the figure drifted considerably 
in the lateral direction as they penetrated the crystal. 

(2) The determination of the range distribution curves for copper 
atoms penetrating the (Oil) surface of f.c.c. copper crystals. Spec- 
ific attention was placed on the effect of impact point on primary 
range and spread. 
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4 . SIMULATION MODEL 



A. Model Discussion 

The original purpose of the GH program was a comparative study of 
binary versus n-body collision processes 0 Their model used a 63 atom 
f.CoC. microcry sta llite . These atoms were capable of movement but 
could not dissipate energy except by further collision „ No surface ef- 
fects were allowed ; which was clearly justified by the short time span 
of the program . An outer shell of immovable atoms was used to simu- 
late a continuation of the lattice and test for complete event contain- 
ment. (This shell is eliminated in the current program), 

The objectives of the GH study were oriented toward analysis of the 
primary atom and one "target” atom. Concentration on the primary’s 
history allowed a streamlining of the program that removed unneeded 
operations and reduced program run time. Essentially^ the computation 
and recording of the motion of the remaining crystal atoms were deleted. 
This deletion was possible since the path of a primary is rarely exposed 
to feedback from target interactions. In analogy , if the path of radia- 
tion damage is represented by an expanding "tree of interactions", then 
the primary is usually spearheading the top branches. 

Another result of the concentrated study on a primary and a single 
target was the reduction in crystal size. It must be made clear that the 
reduction in crystal size has been carefully studied in earlier work by 
Harrison. 15 Considering the uncertainties of the crystal potentials and 
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the use of iteration techniques , the neglect of more distant crystal 
atoms is not significant in regard to primary performance. The crystal 
reduction was carried out experimentally through trial runs on the com- 
puter. The size was reduced until the program run time was optimized, 
but iteraction results still correlated closely. 

The unit crystal used in this study regenerates itself in one of two 
ways. The first form constructs the crystal in the y direction only 
(figure 3) while the second, slower program can construct the crystal in 
any direction (figure 4). As a result of these features, the small unit 
crystal is repeatedly remade in front of the primary. The fast running 
program is used in those cases where the primary is not expected to ex- 
perience severe collisions nor deviate greatly from its original track. 

This type of program is called Chan2. The second program is used in 
areas when the primary does experience large amounts of spread and is 
designated Chan 3. Figures 3 and 4 indicate the unit crystal sizes and 
their orientation. All primaries hit the Oil surface orthogonally (i.e. in 
the y direction). The crystal is positioned by the program with Oil sur- 
face in the x-z plane and with the positive y direction coinciding with 
the primary track at impact. A majority of the runs were made in a tri- 
angular impact area. Several check runs were made in the adjecent triangle 
and a few elsewhere. Detailed discussion of the impact triangles is given 
in the next section. 

o 

The basic building block of a copper crystal is a cube measuring 3.614A 

on a side. The simulation program uses units of length of 1.807A called 
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lattice units . 



The program is written in Fortran and Symbolic Fortran (an assembly 
language closely related to machine language). All runs were conducted 
on a Control Data Corporation 1604 computer which has a storage capacity 
of over 32 , 000 words and six index registers. A CDC 160 computer was 
used to automatically plot graphic projections of the primary trajectories. 

Refer to figure 5 for an example of such tracks projected onto the Oil plane. 

The program constructs the lattice (i.e. assigns positions to each atom), 
commences the run, at designated time intervals computes the forces on the 
primary and secondaries, adjusts their motion accordingly, prints required 
information at set intervals, and "shuts down" when any of three conditions 
exist. The program terminates when (1) the energy of the primary goes be- 
low 25ev, (2) a preset time limit is exceeded (measured in program cycles; 

10 cycles « time for the primary to travel one lattice unit), or (3) the primary 
leaves the lattice (necessary for Chan 2 programs only). The determination 
of this cycle factor received considerable study by Gay^ and was optimized 
experimentally. 

As may have been evident in the brief program description above, the crit- 
icalpha&eof the program is the method of computation of the primary path as 
it experiences crystal forces. In general, the "events" of this simulation 
program are the primary crystal interactions. 

The interaction routine used by GH is still used in this model. It consists 
of a double iteration procedure that calculates the average force on an atom 
over a small distance of travel. This force is used to advance the atom. 
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Harrison states. 

The unbalanced force . . „ (on the primary). . is an average force 
calculated by a double iteration procedure as follows (1) assume 
an atom at position 1 with velocity 1; (2) calculate the total force 
on the atom as a result of all the other atoms in the lattice (this 
means normally only about 8-10 nearest atoms. . .); (3) call this 
calculated force, force 1, and use the equation of motion to move 
the atom to a temporary position, position 2; (4) now repeat the 
force calculations for position 2 , call this force 2; (5) go back to 
positional, and use the average of force 1 and force 2 to move 
the atom to a new position, position 3. Procedures 1 through 5 
constitute one "time step". 14 

The equation of motion used in the iteration procedure differed in two 

1 2 

ways from the equation used by Gibson et. al. Although derived in a 
different manner, the GH equation can be compared more easily to the 
Gibson form by starting with the latter and indicating the modifications. 
The Gibson form of the equation (in the x direction only) is: 

Xi (t + At) = Xi (t) + At[vj (t + A 1/2) + m' 1 . F (t) . AtJ , 
where = initial X coordinate, t = time. At = time increment for change 
of position, Vi = initial velocity, m = mass, and f= the initial force. 
Because At is a small number compared to t the term V* (t + At/2) may be 
expanded in a Maclaurin series. Neglecting all terms after the first , the 
following form is obtained: (t + At ) = xi + At f V i (t) + At/2 . F(2 )/mJ 

In the GH derivation, the force term is an average force where 
- _ Fj + Ff 
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Thus, the GH expression resembles the Gibson expression but utilizes 
a more realistic force term, 

B. Potential Functions 

Two factors dominate all other considerations in the simulation of pri- 
mary interaction in crystals . The first is the general type of model attempt- 
ed. Specifically, does the program use binary or n-body collisions? Is a 
crystal or a random atom arrangement used? The second factor is the se- 
M^iion of a potential function. Many functions have been proposed; none 
have been singled out as best. Until good correlation between experiment 
and computer simulation is made, the proper potential function will remain 
uncertain. In a sense this report and all such simulation studies are an at- 
tempt to locate the correct function for a particular substance. 

Two basic types of potential function were used in this study; Born- 

Mayer and Bohr. The Born- Mayer function used has three parameter varia- 

1 ? 

tions introduced by Gibson et. al. and shown below; 




where r = distance in lattice units, 

a.. 

r o = nearest neighbors distance (2.555A for copper). 



0 



l c = lattice unit (1.804A for copper). 



p = input parameter, and 



A = input parameter (ev). 
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The form shown above is convenient for computer calculation. The 
three sets of parameters selected by Gibson are? 



Potential 


A(ev) 


0 


I 


.0392 


16.97 


II 


.051 


13.00 


III 


.1004 


10.34 



To speed calculation , the Born-Mayer force functions were truncated at f=10 

o 

newtons (corresponds to approximately 3A) . 

In similar format the Bohr potential function is: 

0b = [Vr] exp (On c b a b - r/a b ), 
where r = distance in lattice units, 

a b = input parameter in lattice units, and 
c b = input parameter in ev. 
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The two parameters are % c b = 9.94x10 ev, and a b = .06741u. These para- 
meters are the ones used in the RO report. Further, the corresponding force 
truncation value was used so that the Bohr potential functions and force func- 
tions in this study are identical to Robinson's. (The truncation of the force 
function occurs at f = 3.9x10“° newtons). 

The lattice has zero kinetic and potential energy prior to primary impact. 
In other words no vibrational energy is present. 

C. Binary vs n-body model 

This study made use of an n-body simulation program in contrast to the 

113 17 

simplified binary collision approach of Robinson, ' ' et. al. The in- 
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herent danger of superimposing independent collisions is that simultaneous 
collisions will not be properly accounted for 0 For purposes of fast comput- 
er runs, the binary program has the greatest advantage, but its justifica- 
tion requires better experimental verification than is yet available. The 
n-body program is inherently more accurate because of its more thorough in- 
teraction routine. Given proper parameters, it should yield the more accurate 
data. 
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5 „ EXPERIMENT PROCEDURE 



Computer simulation studies are experiments in the accepted sense of 
the word. They are trials made to confirm or disprove some suggested 
phenomenon, and are subject to most of the usual experimental inaccuracies „ 
Systematic errors occur because of improper model simulation or necessary 
oversimplifications of the model. Random errors caused by computer mal- 
function should be rare. Provided the program and its input data are not 
modified, the computer should reproduce results perfectly. (Refer to subsec- 
tion C.). 

Simulation studies are sensitive to parameter variation to the extent that 
some computations can become meaningless if not properly controlled. For 
example, the combination of high energy primaries, Bohr potential, and sel- 
ected impact points produce an interaction so weak that the program cannot 
maintain the energy check. As a result the primary could penetrate infinite- 
ly. Such extreme situations have been avoided and are not included in the 
data of the report, but they do point out the danger of blind faith in computer 
"answers" . 

The remainder of this section will present the experimental approach used. 
A. Determination of impact area and sample size 

Symmetry considerations reveal that the impact rectangle shown in figure 
3 or 4 represents all possible combinations of impact possibilities. That is, 
the rectangle is a "unit area" that, if repeated, will completely describe the 
Oil surface. Further, with one exception, the rectangle can be split down 
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the diagonal into two triangles that have point correspondence in regard to 
primary interaction. Consider impacts at points A' and C'(the 90 degree an- 
gles). As indicated in figure 4, atom row AB is nearer the surface than row 
DC and therefore "sees" the primary sooner. Now, a primary, or "bullet", 
striking the crystal at A' "sees" these top plane or near atoms (row AB) at 

o 

a distance of 1.28A . However, a bullet hitting at point C' sees the top 

o 

plane atoms (row AB) at a distance of 1.807A . Obviously, this difference 
effects the first interaction and subsequently alters the bullet” s path and 
future collisions. It will now be shown that this difference in the two 
triangles, will cause a mild perturbation in range and spread data, and is 
not considered serious for purposes of this study. (The triangular impact 
area is desirable since it reduces the sample size necessary by one half). 

A graphical correlation analysis of two critical parameters (range and 
spread) was performed to insure the validity of the impact triangle. In 
each triangle 14 points were chosen and the range and maximum spread 
values recorded. All points were obtained using lOOOev primaries and a 
Gibson II potential. Since this combination has a strong primary-crystal 
interaction, discrepancies between the two triangles should be maximum 
(figure 7). In general, good correlation was revealed. Note that the range 
correlation values show less correlation for deeper penetrations. This is 
expected since those primaries that impact further from the "corner atoms" 
have smaller energy loss rates and more chance to deviate. The spread 
terms have less correlation at small values. However, these terms are 
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sufficiently small so that .2 or .3 lattice unit variation between triangles 
does not effect the spread contour diagrams (to be discussed later). As a 
result of these considerations, triangle ABD was selected as the primary 
impact area. The error introduced by this assumption is small and is not 
known precisely. In future studies one area of further investigation should 
undoubtedly be a more thorough check of the rectangular impact area and its 
effect on results . 

From the grid shown in figure 8, 14 impact points were selected and 
are labeled one through 14 . For some of the runs these original 14 points 
were supplemented by 12 more. These points are labeled in the figure as 
1A through 12A. The supplemental 12 point array was chosen to verify the 
selection of the basic 14 point sample size . Both the Gibson II- 1000 and 
5000ev studies were initially made on 14 point grids and later supplemented 
by 12 more data points. The range contour surfaces remained regular when 
the new points were added, and showed no anomalies (figure 9). Further, 
the range distribution curves for the 14 point and 26 point grids correlated 
very closely (figure 6), and established the soundness of 14 point grids. 
These range curves were constructed on an area rather than a point basis. 

That is, cross sectional areas of the contour surfaces (figures 9 through 12) 
were measured at specified depths of penetration. The ratio of this area 
to the total triangle area gives the percent of primaries stopped. In figure 
9 the shaded area represents those primaries stopped within 100 lattice units. 
The ratio of this area to the total area provides the percent stopped at 100 
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lattice units. (The range curves were constructed from larger contour 
diagrams than those shown). The area method yields more accurate pene- 
tration curves than would a point system. Overall justification of a 14 
point sample is simply the general regularity of the results. 

B. Impact points A,B,C and D« 

In addition to the 26 point grid, four other points were chosen on the 
triangle. These points coincide with the impact points in figure 3 of the 
RO study (three of the four points showed spread). These points are 
shown in figure 8 of this report and are labeled as A,B,C and D. 

As an initial step in this study, the parameters of the RO "A,B,C 
and D shots" were duplicated as closely as possible. Later, several 
other combinations of bullet energy and potential were used. Table 1 
shows the number of grid points investigated under each combination of 
parameters. As will be shown in the next subsection, many of these 
points were rechecked with additional runs. 
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TABLE 1 





Impact Points Studied By Parameters 








RO 13 


Points 


14 Point 


12 


Point 


Symmetrical 




A,B 


,C,D 


Grid 


Grid 


Triangle 




Energy 


Energy 


Energy 


Energy 


Potential 


1000 


5000 


1000 


5000 


1000 


5000 


1000 


Bohr 


4 


Point C 


14 


13 


8 


6 


- 


Gib II 


4 


Point C 


14 


14 


12 


12 


14 


Gib III 


4 


- 


14 


- 




- 


- 



C. Reproduced bility of runs. 

As stated earlier reproduced bility problems should not occur in com- 
puter studies. Unfortunately, the program - computer combination used 
in this study does fail to reproduce periodically. As a result a large 
number of reruns were necessary. The sample was random, but was 
later checked to insure that most combinations of energy, potential and 
run length were represented. The sample size was approximately 15 per- 
cent. It was found that the probability that a run would be reproduced 
exactly was 85 percent. Further, once a run was confirmed by a rerun, 
the probability cf its correctness could be assumed to be one. The lack of 
perfect reproduceability has not been determined. However, it was found 
that those runs that deviated seemed to get out of phase, but did not vary 
greatly in amplitude or distance traveled. Generally, the variation of am- 
plitude (i.e. spread) and depth of penetration were less than a few per- 
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cent. In one case the depth of penetration varied by about 5 0 percent, 
but the depth of penetration in this case was very short ( <8Lu). 

D. Program output. 

The output will be briefly discussed here, so that subsequent dis- 
cussions of data will have more meaning. The program as used in this 
study was designed specifically to examine only the primaries' per- 
formance. Output is made every ten cycles. Ten cycles is the time 
necessary for the bullet to travel approximately one lattice unit (1.807A). 
Each output yields the following data about the primary; 

(1) x,y, and z distances (measured in z lattice units) from the im- 
pact point (a "z lattice unit" is 1.414 times greater than either the x or 
y lattice units) , 

(2) x,y, and z velocity values in "z lattice units" per second, 

(3) kinetic energy, potential energy, and total energy in electron 
volts, and 

(4) real elapsed time in seconds. 

Except for a few 6000 cycle runs all of the simulations used a 2000 
cycle cut off. 
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6. RESULTS AND CONCLUSIONS 



The data collected from the printout sheets discussed in section 
five were condensed and tabulated as follows: (1) primary final energy 
(either approximately 25ev or the energy at the 2 000 or 6000 cycle cut 
off), (2) depth of penetration in lattice units, and (3) the maximum 
spread. Recall that spread is the perpendicular radial distance (in 
lattice units) that a primary deviates from its initial straight line tra- 
jectory. 

The type of tracks reported by RO (those of large penetration and 
spread) suggest the need for a word descriptive of their nature. There- 
fore, in this report, the term wander will apply to all primaries that ex- 
ceed two lattice units of spread, but whose energy loss rate does not 

o 

exceed 25ev/A. 

Four combinations of potential and energy were studied: Bohr and 
lOOOev primaries, Bohr and 5 000ev primaries, Gibson II and lOOOev 
primaries, and Gibson II and 5 000ev primaries. Since the data was 
collected from the impact point grids, it is best displayed using a pic- 
torial representation of the grid. The penetration values are shown at 
their respective grid impact points and can be visualized as a penetra- 
tion contour surface (figures 9 through 12). The pertinent data from the 
spread contours was summarized and is discussed below. 

A. Spread results 

The spread values for each impact triangle are a function of primary 
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energy, potential and impact point. Separate spread contours are in- 
cluded in the report, and are located on the same page as their corres- 
ponding penetration contours. 

The results of the spread data collected are summarized below? 

(1) Wander, as observed in the Bohr lOOOev primary case by RO, was 
found to exist for several other combinations of primary energy and poten- 
tial function. 



TABLE 2 



Percent of primaries that exceeded 2LU spread 
and Percent of primaries that wandered 




POTENTIAL 

Gib II Bohr 


Energy 

lOOOev 5000ev 


Energy 

lOOOev 5000 ev 


Percent 
>2LU Spread 


27 


54 


43 


25 


Percent 

Wander 


0 


11.5 


5 


15 



The result of RO can be generalized? When parameters are such that 
penetrations vary from deep to short values over the surface of impact 
triangle, then some of the primaries will wander. Note that the Gibson 
II - lOOOev series of runs does not satisy this criteria (refer to table 2), 
since all the primaries had short penetrations. 

From another viewpoint, if the strength of interaction between primary 
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and crystal varies over the impact grid from hard to soft collision, then 
some primaries will wander. Thus, the presence of wander is not simply 
a result of the Bohr potential, but is a function of the strength of inter- 
action. 

(2) A study of the spread data (figure 9 through 12) reveals that the 
primaries that wander are located in a particular region of the impact 
triangle. Note that the primaries impacting in zone A, figure 14, have 
large spread values while those in zone C have large penetrations. It 
is not surprising that the primaries that wandered were located at the 
intersection of zones A and C. The overlap of the two zones varies some- 
what with the parameters of the run and cannot be firmly located . How- 
ever, the important result is that wander is a function of impact position 
and is simply the transition region from the short range penetration zone 
to the deep penetration impact zone. 

(3) The outlined rectangular area in zone C (figure 14) has primaries 
that have low spread. Any atom impacting this area will have large pene- 
tration and will not exceed .8 Lu spread. The primaries in the remaining 
area of the impact triangle have at least lLu of spread except for a 
small area near point 11; the midpoint between the atom crystals. This 
midpoint is a symmetrical force region for the primaries and little deflec- 
tion or spread takes place. Further, in the Bohr potential cases, the force 
curve is near cut off so that the forces on the primary are negligible. Two 
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regions of large spread occur between the ,3 and . 5Lu radii of the cor- 



ner atoms. These are the zones where the primaries wander. 

(4) Primaries 2, 5A, and 8A in figure 15 wandered. As shown by 
the "history curves" of the runs, the rate at which energy is lost is us- 
ually uniform until approximately 1000-1 2 OOev and then the energy loss 

is more rapid. The general appearance of the family of curves is uniform, 
and one of gradual transition from short to long penetration as the impact 
points move from the corner atoms of the triangle to the channel area 
near the right angle. The short penetration curves display less regularity. 

o o 

The gap in the region of 180A to 300A penetration is not an amonaly. 

If primaries were concentrated in the region of rapid penetration change 
(i.e. the "hills" of the penetration contour) then this gap would be fill- 
ed. Though the crystal effect on range and spread seems logical, the 
rigid consistency and smooth transition is surprising. For example, if 
more impact points are studied, the occurance of occasional penetrations 

o o 

in the range of 300A to 600A that do not wander can be expected. 

In brief summary, primaries do wander if the impact point is located 
in a transition region from hard interaction to soft interaction. The "wan- 
der phenomenon" is part of a smooth change from one type of interaction 
to another. 

(5) The four impact points of the RO study that wandered were investi- 
gated (figure 2). Their parameters were duplicated (Bohr truncated poten- 
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tial, lOOOev) to insure a valid comparison 0 In addition, the Gibson II 
and Gibson III potentials were used. The selection of the exact para- 
meters in the Bohr case provided a comparison between the RO binary 
model and the n-body model. 

Of the four points, only point C displayed wander and it did so with 
the following combinations of paramenterss (1) Bohr and lOOOev pri- 
maries, (2) Bohr and 5000ev primaries, and (3) Gibson and SOOOev pri- 
maries. As seen in figure 14, points A, B, and D are in zone C and 
are typical of that zone since they displayed deep penetration and less 
than ,8Lu spread. The probable reason for the disagreement between 
the RO study and this study is the use of the n-body model. (Refer to 
the next subsection for detailed comparison of the results obtained from 
the two models .) 

At present, none of these results can be verified experimentally. If 
the potential functions studied are not reasonable, the conclusions dis- 
cussed above could be radically change. As will be seen in the next sec- 
tion, the potential functions selected (particularly the Bohr) may be too 
weak. If so, the amount of "wander" found in future studies using "hard- 
er" potentials may decrease or disappear. 

B. Penetration results 

The penetration results are shown in contour diagrams (figure 9 through 
12) and the integral penetration curves of figure 6. Penetration contours 
for all runs are regular and well defined. The following points are reveal- 
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ed from the data t 



(1) Impact points 4, 4A ,7 and 7A, which are located in the 110 
channel of the crystal have maximum penetration „ The Gibson II- 
lOOOev contour is complete and reveals the smallest and lowest 

peak of maximum penetration,, The Gibson II-5000ev contour is virtual- 
ly complete (four runs truncated at 600Lu) and is shaped similarly to 
the Gibson II lOOOev contours 0 However, the contour height is much 
greater and the region of very low penetration (below 20Lu) is much 
smaller. 

Both Bohr contours were truncated at 2000 cycles „ In comparison to 
the Gibson results, the areas of low penetration are smaller and the rate 
of change to deeper penetration more rapid. Excessive computer run 
time curtailed the construction of complete contour surfaces for the Bohr 
potentials . 

The Bohr- lOOOev contour and the Gibson II~5000ev contour are simi- 
lar in size and shape, which indicates that these combinations of ener- 
gy and potential produce similar interactions. 

(2) As the energy of the primary is increased the depth of penetra- 
tion increases; also, the amount of penetration varies with the potential 
function used. 

The n-body integral penetration curve (5000ev and Gibson II poten- 
tial has much larger penetrations than the RO results (figure 6). Un- 
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fortunately, none of these simulation curves can be easily compared 
to the experimental data of figure 1, because of the large difference 
in the substances studied. Important conclusions are available from 
the results, however. 

The 5000ev curves in figure 6 of RO and this study agree fairly 

o 

close in the region of small penetration (up to 120A). Thus, it is pro- 
able that the n-body program agrees closely with the binary model for 

o 

severe interactions (i.e. small impact parameters below . 3Lu or . 5A ) . 

A thorough study by GH established the close agreement between single 
binary collisions and n-body collisions at short impact parameters. 

Their study investigated a wide range of energies and several poten- 
tials (including the Gibson II). 

(3) In the region of deep penetration the two models disagree . 
Nevertheless, an important fact was revealed when the simulation pene- 
tration data were compared with the corresponding potential curves. 

Refer to figure 13 which shows the Gibson II potential discussed here, 
the Bohr potential, and a new potential to be explained in the next sub- 
section. As seen in the diagram, the three curves in question use the 
basic Gibson II potential function. However, RO curves B and C are 
truncated at much larger potential values than the potential of this study. 
As seen from the integral penetration curve and these potentials, as 
the distance to the truncation point is increased, the depth of penetra- 
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tion of the primaries increase . (This effect is opposite of what might 
be first supposed). The most plausible explanation probably lies in 
the subtle effect of the small potential curve slope (hence small force 
values) that occured in the region beyond the truncation of the RO 
curves. These small forces confine the primaries to the more open 
crystal regions. As an analogy, the deeper penetrations can be com- 
pared to a rock skipping over the surface of water. It makes many 
weak impacts that confine its movement. This is the case of the n- 
body potential. The more radically truncated potentials of the RO study 
probably cause abrupt interaction when a primary moves from a region 
of no potential to one of large finite value. Thus, deep penetration of 
primaries still occurs in the RO model but to a much smaller extent. 

In summary, one significant factor causing the deeper primary pene- 
trations of this study was the use of a potential that was truncated at 
large radii. The use of the n-body model also contributed to large pene- 
trations and will be discussed shortly. 

(4) As shown in figure 13, the RO study also used an eroded poten- 

o 

tial that decreased to zero at 1.27A. The corresponding integral pene- 
tration curve is similar in appearance to the n-body penetration curve of 
this study, but has slightly lower penetration values. Thus, the poten- 
tial function if eroded or truncated at large distances tends to produce 
similar large penetration depths. In contrast, potentials that have signifi- 
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cant values at truncation reduce penetration depth. 



(5) Another factor, which undoubtedly caused the deeper penetra- 
tions found in this study, is the n-body model* It s use insures a 
continual superposition of all significant crystal forces; and consequent- 
ly, the net force on the primary is usually smaller and produces less 
scatter. Thus, the wander detected by RO for impact primaries A,B, and 
D, but which was not verified using the n-body model, may result from 
the insufficient superposition of forces in the binary approximation collis- 
ion process. This helps clearify the results in the previous subsection 
concerning wander. 

C. Recommendations 

Two proposed simulation studies are presented in this subsection. 

Both are attempts to obtain verification with experimental data. Recent 
18 

studies by Abrahamson provide potential functions for tungsten and the 
inert gases. This development, in conjunction with the experimental 
work of Davies, provides, for the first time, the possibility of positive 
experimental verification of computer simulation programs. 

It is recommended that the n-body program incorporate the new poten- 
tial function developed by Abrahamson, and that it be used to simulate 

135 

the Davies experiment of Xe in tungsten. 

Although the experimental results for copper primaries in copper are 
not available, it is felt that the simulation penetration values are too 
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large (especially the n-body) and should agree more closely with the 

25 kev - Kr 84 in Cu experiments and the Xe 135 in Cu experiments at 

5kev. In this regard RO have used "Cu" primaries of triple mass 
192 

(Cu ) which did reduce penetration ranges (fig. 6), but not suffi- 

84 

ciently to agree with the results of 25kev Kr in Cu„ 

In spite of the probability that the n-body integral penetration curve, 
obtained in this study, is further from experimental fact than the RO 
curve, the n-body model is an inherently more accurate approach. This 
conclusion is based on previous discussion which reduces to the essen- 
tial fact that the n-body model is more realistic and, hence, a better 
simulator „ Equally critical , the results of the RO study were obtained 
using severe truncation values (fig 13). Truncation of the force func- 
tion at significant values is unrealistic and probably leads to incorrect 
results. To obtain closer verification of experimental results, stronger 
potential functions may be necessary in future n-body program studies. 
(It is possible that experimental integral penetration curves could be 
verified by most models by sufficient adjustment of the potential 
functions. Thus, experimental verification in a single area will probab- 
ly not prove sufficient, in the long run, to establish a particular poten- 
tial function ) . 

Preliminary calculations indicate that the potential function of 
Abrahamson (Xe in W) is significantly stronger than the Gibson II poten- 
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tial. If so, it is likely that this potential will produce excessive 
duction in primary penetrations and its use in a simmulation program will 
not yield close correlation to Davies experimental work 0 

It is also recommended that the verification of Davies experiment 
be attempted using a "compound" force function. This function should 
follow the Gibson II force curve at small impact parameters where 

o 

good agreement is obtained at present. At approximately 1. OA the 
curve would change to a different exponential form that exerts much 
larger long range forces on the primaries. This function will reduce 
deep penetrations and agree more closely with experimental data. The 
parameters of Davies Xe in W should be duplicated for this study also. 



34 



Res idual Activity 



Key 



— Xe^'" in w (Davies et* al.^) 

— Kr^’ in Cu (Lutz et. al. 10 ) 
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Fig. 2. ROBINSON AND OEN ,r . ..(Oil) channel trajectories onto (Oil) i3 
surface for f.c.c. Cu. ..(lKev Ca primaries ) . . . Bohr potential."' 
Points A,B,C and D penetrated 300A in part of trajectory shown. 
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Near atoms 









2 



(Z)^- 



Far atoms 




1 2 3 

X( Lattice Units) 

Fig. 3 . (HO) FACE OF ’ CHAN2’-. TYPE CRYSTAL. Shaded triangle indicates 
impact area; (*) RO^ ’'wander' 1 impact point coordinates. 
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1 2 3 U 5 6 

X (lattice Units) 

Fig. li. (110) FACE OF CHANlt TYPE UNIT CRYSTAL. Shaded triangle indi- 
cates impact area # 
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B 
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D 

E 

F 



Model 


Primary 
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N-body 

Binary' 1 -* 


5 

5 


Binary^ e 


5 (amu=6U) 


Binary^ 


5 (SKU- 193 ) 


N-body 


l 


N-b’ody 


l 



Potential 
Function 
Gibson H 
Born-Mayer II 
Born-Mayer I 
Gibson II 
Gibson H 
(lii pt. sample) 
Gibson II 
(26 pt. sample) 



Fig. 6. INTEGRAL PENETRATION CURVES. 

(Curves E and F are essentially superimposed.) 



4 0, 




Symmetrical Scatter (Lu) Symmetrical Penetration (Lu) 




Impact A Penetration (Lu) 




Impact 'A Scatter (Lu) 

Fig. 7. DEPTH OF PENETRATION AND MAXIMUM SCATTER CORRELATION 
PLOTS (impact A vrs. symmetrical A). 
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rcA- 3 "wander" impact points. 





Fig. 9. CONTOUR CURVES. Gibson II potential and 5>000ev primaries. 
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Maximum Scatter Contour (in Lattice Units) 




Penetration Contour (in Lattice Units) 



Fig. 10. CONTOUR CURVES. Gibson II potential and lOOOev primaries. 
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Maximum Scatter ,Contour (in Lattice Units) 





Penetration. Contour (in Lattice Units) 

Fig. 11. CONTOUR CURVES* Bohr potential and £OOOev primaries. 
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Penetration Contour (in Lattice Units) 



Fig. 12. CONTOUR CURVES. Bohr potential and lOOOev primaries. 
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Zone Ar Region. of large primary spread. 

Shaded area describes region in which wander may occur. 

Zone’ C: Region of large primary penetration. (Zones do not apply to lOOOev 
Gibson II data.) 




(A0) jQaaua jCjeuipaj 
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Fig, 15, GIBSON II POTENTIAL, 5000ev primaries. Indicates rate of loss of 
energjr of primary as function of impadt point. 
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APPENDIX II 



PROGRAM DISCUSSION 

A. Program Explanation 

The descriptions given below are to be used with Appendix I (Pro- 
gram Listing); and are specifically applicable to programs Chan3, Bull- 
ets , or G Chan3 . 

BOX 1 

Listing of dimension # common; and Format statements. Provide graph 
plot data needed such as graph size, scale, and titles. Read force and 
potential function parameters. Read primary atom parameters. Evaluate 
miscellaneous constants. 

BOX 2 

Designation of nearest neighbor lattice distances. Construction of 
microcrystallite. Three nested "DO loops" that extend to statement 63 
are used to locate and number the atoms in the crystal. 

BOX 3 

Initialize the velocity and the position of the primary. 

BOX 4 

Store initial positions of all atoms. Initialize various terms (i.e. 
set=0). Print variable quantities such as primary energy, direction 
cosines, mass, etc. which comprise the output heading. 

BOX 5 

Compute several constants given DT, CVD, PTMAS and PGMAS. 
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BOX 6 



This entire box is a "DO ioop" on all crystal atoms. Set all forces 
to 0, Check coordinate distances from primary to lattice atoms and 
neglect if greater than ROE, Compute the square of the "vector" dis- 
tances between the primary and lattice atoms and neglect if greater than 
2 

(ROE) . Compute the force that corresponds to the vector distance and 
compare it to FM. Neglect if less than FM. Compute forces on the 
primary and lattice atoms , 

BOX 7 

If index=l jump to box 8. Compute temporary position of the primary 
and then the rest of the crystal atoms. All atoms are advanced at their 
initial velocity for DT seconds, while experiencing forces found in box 6. 
BOX 8 

(Return to box 6 and recompute the forces on all atoms. These values 
are added to the initial forces yielding force terms twice the average. 

Since index=l on reaching box 7 a jump is made to statement 260 in box 8). 
The averages of the initial and temporary position forces are used to compute 
the final positions of the primary and lattice atoms. (See page 14 of text). 
BOX 9 

Perform tests to see if program output is necessary or if the program 
should rebuild the crystal. If the latter is true jump to the appropriate re- 
generation box. 
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BOX 10 



Calculate the distance from the primary to each lattice atom; if less 
than ROE calculate the potential energy of the interaction and store. 

BOX 11 

Calculate the quantities needed in the output statement (i.e. pri- 
mary position, kinetic energy, total energy, and real elapsed time). 

Write the output statement. Check to insure the potential energy is less 
than "Bitty". Terminate program if the potential energy is less than "Bitty" 
for three passes. 

BOX 12 

Check current program cycle. If at "cut off" cycle, terminate the 
program. 

BOX 13 

This section regenerates the microcrystallite in the+Y direction. 

BOX 14 A 

Regenerate microcrystallite in the -X direction. 

BOX 14B 

Regenerate microcrystallite in the +X direction. 

BOX 15A 

Regenerate microcrystallite in the -Z direction. 

BOX 15B 

Regenerates microcrystallite in the +Z direction. 
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BOX 16 



Call the graph plot subroutine (called Draw ) . The following parameters 
are necessary: 



Variables 


Explanation 


II 


Number of points to be plotted 


GRX 


X coordinate 


GRZ 


Y coordinate (simulation program 
Z coordinate) 


C 


Number of curves per graph 


0 


Point or curve plot 


LABEL 


Type End by the end of curve 


ITITLE 


Graph title 


EXSCALE 


Graph scale in X direction 


YSCALE 


Graph scale in Y direction 


0 


Distance of x axis from bottom 
of graph 


0 


Distance of y axis from left edge 


2 


Mode of x axis (automatic) 


2 


Mode of y axis (automatic) 


9 


Graph width 


15 


Graph heighth 


1 


Place grid on graph 


ERROR 


Indicates if previous plot was OK 
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B. Definition of Variables: 



Term 


Usaqe 


Bitty 


a small increment that is compared to 
primary potential energy (PTE) . Used 
to terminate program when potential 
energy goes below increment value. 


BE 


primary energy (in joules) 


BX, BY, BZ 


primary initial coordinate (in lattice 
units) 


COX, COY, COZ 


direction cosines of primary 

o 


CVD=1.8075xl0' 10 


converts lattice units to A (Cu only) 


CVE=1.60xl0 -19 


converts ev to joules 


CVM=1 . 672xl0~ 27 


converts AMU to kg. 


DIST 


distance between two atoms or the 
square of the distance (lattice units) 


DT 


basic time step length (seconds) 


DTI 


variable scale factor for DT 


DTOD 


DT/CVD (a "units converter") 


DRX, RRY, DRZ 


coordinate d istance of primary to 
lattice atoms (in lattice units) 


DXS , DYS, DZS 


locates current microcrystallite with 
respect to original origin (lattice units) 


DXT, DYT, DZT 


test constants used to rebuild crystallite 


EI=25 . 0 


cut off energy (in electron volts) 


EXA, EXB, EXC, EX1 


input potential function parameters 
(see Gay 1 ^ thesis) 


EV 


primary energy (in electron volts) 
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FA t FB , FC 


x,y, and z force components be- 
tween two atoms (in newtons) 


FM 


a small force increment used to eli- 
minate weak force terms 


FOD 


force/distance (in newtons/lattice units) 


FORCE 


eroded force between two atoms 
(in newtons) 


FXA , FXC 


input force function parameters 
(see Gay thesis 


FX, FY, FZ 


x,y,z coordinate force terms, respectively 
(in newtons) 


GMAS 


primary mass (in AMU) 


GRX, GRZ 


coordinates of primary that are used by 
graph plot routine (in lattice units) 


IDX, IDY, IDZ 


constants used in lattice reconstruction 


IT, JT, KT 


test points in lattice generator 


DC, IY, IZ 


unite lattice dimensions (in lattice units) 


IXP, IYP , IZP 


specify volume of microcrystallite (in 
lattice units) 


IYT 


used in crystal regeneration 


HDTOM 


DT/ (2*PTMAS) 


HDTOMB 


DT/(2*PGMAS) 


HGMAS 


!/2 primary mass, units: (ev*kg)/ 
(m 2 / sec 2 ) 


HTMAS 


1/2 target mass, units (ev*kg)/ 
(m 2 / sec 2 ) 


M=2 


a constant 
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MASK 


used to calculate absolute value 
(removes minus signs) 


ND 


print statement increment (in cycles) 


NS 


initial print statement (in cycles) 


NT 


current program cycle 


NTT 


program cycle for final print statement 
("cuts off program") 


ONE=l . 0 


a constant 


OPNX 


1 . 9 SOX 


OPNZ 


1.9 SCZ 


PGMAS 


primary mass (in kg) 


PKE 


primary kinetic energy (in electron volts) 


POT 


potential between any two atoms (in elec- 
tron volts) 


PPE 


primary potential energy (in electron 
volts) 


PTE 


primary total energy (in electron volts) 


PTMAS 


target mass (in kg) 


QDTOM 


DT/ (4* PTMAS) 


QDTOMB 


DT/(4*PGMAS) 


ROE 


nearest neighbors separation (in lattice 
units) 


RO 


= 1/ROE 


RX, RY, RZ 


atom coordinate (in lattice units) 


RXI, RYI, RZI 


initial atom coordinates (in lattice units) 
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RXK, RYK, RZK 


temporary atom coordinates used in 
force calculation (in lattice units) 


SCX, SCY, SCZ 


factors to convert X,Y,Z to the same 
lattice unit basis; (values depend on 
type crystal used; i.e, BCC or f.c.c. 
etc. ) 


TIME 


real elapsed time (in seconds) 


TMAS 


target mass in AMU 


TPO=2.1 


a constant 


TPOX 


2 . 1 SCX 


TPOZ 


2.1 SCZ 


VOL 


velocity of primary (meters/seconds) 


VX, VY, VZ 


atom velocity (in meters/seconds) 


XLL, ZLL 


test constants that determine if 
regeneration is necessary 


ZE= 0.0 


a constant 
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APPENDIX III 



PROGRAM OPERATING PROCEDURES 

The simulation model consists of a single basic porgram that has 
two forms: 

(1) A fast running program that can reconstruct the basic lattice in 
the y direction only (the direction of the incoming bullet) . This pro- 
gram is limited to cases in which the bullet does not intereact strong- 
ly and hence doesn't leave the basic lattice in any direction except 
the y direction. This program runs in about one half the time of the 
other program and has three identical decks which are titled as follows: 
Chan, Gchan, and Bullet, For runs cut off as 2000 cycles the "Chan 
type " program requires about ten minutes. Those shots that strike 

the impact triangle near the corner atoms often terminate in very short 
time since the primary rapidly drops below the cutoff energy of 25ev. 

(2) The slower program can contain any bullet trajectory since it 
reconstructs the lattice in any direction required. The testing procedure 
to perform this function increases the program run time. This program 
has two identical decks that are titled: Chan3 and Bmlet3. (The various 
deck names serve only to expedite processing and ease administration 
of runs. The slower program (Chan3) was the third modification of the 
basic n-body program). Run times are about 20 minutes for 2000 cycle 
cases . 

The items discussed in the remainder of this section are applicable 
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to all decks : 



A. Parameter cards internal to the deck 0 

Within the program deck there are three parameter cards that are often 
changed; the potential function, force function, and the force truncation 
card. The force and potential cards are the third and fourth cards in all 
decks . These fortran functions are computed from a single imput data 
card (see the next sub section). Thus the force and potential cards are 
paired to an input parameter card. Specifically, the Bohr potential and 
force cards have a single corresponding input card. The Born- Mayer 
force and potential cards become one of the three Gibson forms (1,11, or 
III) through the use of one of three corresponding input parameter cards. 
The third internal card is the force card (ex. fm = 1.0E-11) which allows 
the selection of any desired truncation force. 

B. Input data cards » 

These are the basic data cards used at the end of the deck to spec- 
ify the various parameters for a run and to set up the graph plot routine. 
The following list describes the cards in order as they would appear for 
a single run: 

(1) The first card prints any desired statement at the top of the 
print out. 

(2) This is the force and potential input parameter card and it supplies 
the follwoing values: EX1 , EXB , EXC, FXC, IH. The first four quantities 
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are needed to define the function and the last (IH)prints the potential 
function in the printout „ The above data cards are used once per run. 
(Chan5 uses a new input parameter system )„ 

(3) This card provides the variable quantities that specify the 
primary parameters and general program information „ The following 
items appear: 

GMAS = Primary mass in AMU 
TMAS = Lattice atom mass in AMU 

BX = X initial coordinate for the primary 
BY = Y initial coordinate for the primary 
BZ = Z initial coordiante for the primary 
EV = Primary energy in ev 
COX = Direction cosines of primary 
COY = From the Oil plane to trajectory line 
NTT = Cut off cycle (ends program) 

NS = First print statement (in cycles) 

ND = Print increment (in. cycles) 

Recall that one cycle = . 1 time for the primary to travel one lattice unit 
(1.807A). 

The next six data cards are used to title and identify the graph out- 
put for the run. 
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(4) Potential is designated (ex. GIB. for Gibson) 

(5) Energy (ex. 5000ev) 

(6) X coordinate (ex. x = 2. 00) 

(7) Z coordinate (ex. z = 3.00) 

(8) Cosine X (ex. cox = .00) 

(9) Cosine Y (ex„ coy = 1.00) 

This completes the data cards needed for single run. 

To provide additional runs cards three through nine are repeated for 
each run. An "end card" is placed last on the deck to terminate the pro- 
gram. The program is designed to run on the regular monitor operation 
service, but for runs longer than 15 minutes (the operation service limit) 
the program must be run individually. Three tapes are required, in add- 
ition to the Fortran 60 compiler tape, to run the program on the 1604 
computer. These tapes are the input, output, and graph plot tapes. 

(To perform runs without the graph plotter pull the "call graph card", 
the title cards (4 thru 9), and the card that "reads in" the title cird. 
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